Theoretical analysis for critical fluctuations of relaxation trajectory near a 

saddle-node bifurcation 
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A Langevin equation whose deterministic part undergoes a saddle-node bifurcation is investigated 
theoretically. It is found that statistical properties of relaxation trajectories in this system exhibit 
divergent behaviors near a saddle-node bifurcation point in the weak-noise limit, while the final 
value of the deterministic solution changes discontinuously at the point. A systematic formulation 
for analyzing a path probability measure is constructed on the basis of a singular perturbation 
method. In this formulation, the critical nature turns out to originate from the neutrality of exiting 
time from a saddle-point. The theoretical calculation explains results of numerical simulations. 
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I. INTRODUCTION 

To uncover the nature of fluctuations near a bifurcation 
point has provided a clue to understanding of singulari- 
ties observed in a rich variety of phenomena. The most 
typical example of such studies is the Ginzburg-Landau 
theory for equilibrium critical phenomena Accord- 
ing to this theory, the description of fluctuations in the 
system that undergoes a pitchfork bifurcation is a start- 
ing point for characterizing the Ising universality class of 
the paramagnetic-ferromagnetic transition. The second 
example is a theory of collective synchronization in cou- 
pled oscillators @. In this phenomenon, the mechanism 
of the cooperative oscillation is explored by studying fluc- 
tuations near a Hopf bifurcation. The third example is 
related to a theory of directed percolation, whose univer- 
sality class is characterized by a transcritical bifurcation 
with a multiplicative noise Q . From a viewpoint of bifur- 
cation theory, a pitchfork bifurcation, a Hopf bifurcation, 
and a transcritical bifurcation are local co-dimension one 
bifurcations 0| ■ Now, the last one in this type of bifur- 
cations is a saddle-node bifurcation. 

Mathematically, a saddle-node bifurcation in a differ- 
ential equation is defined as the appearance of a pair of 
saddle type fixed point and node type fixed point with 
respect to change in a system parameter. This bifurca- 
tion has been found in many models such as a mean field 
model of the spinodal transition [5], a model of driven 
colloidal particles @, bio-chemical network models @- 
[To| . a dynamical model associated with a /c-core perco- 
lation problem [TTJ] , and a random- field Ising model [l2j . 
Theoretically, once it is found that a system undergoes a 
saddle-node bifurcation, its deterministic behavior near 
the bifurcation point is immediately derived, as seen in 
standard textbooks [4|. As an example, a characteris- 
tic time scale exhibits divergent behavior proportional to 
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e -1 / 2 , where e represents the distance from the bifurca- 
tion point. 

Now, in a manner similar to the other local co- 
dimension one bifurcations, it is expected that fluctua- 
tions near a saddle-node bifurcation exhibits a singular 
behavior. Thus far, the singular nature of the fluctu- 
ations has not been focused on except for some works 
[E H, S ITll4l3 | , and its theoretical study is still imma- 
ture compared with much progress in theories of fluc- 
tuations near the other bifurcations. However, as we 
already pointed out [H, Illl - ll3j |. a stochastic model un- 
der a saddle-node bifurcation is regarded as the simplest 
one of systems that exhibit the coexistence of discontin- 
uous transition and critical fluctuation. Such a coexis- 
tence has been observed in the dynamical heterogene- 
ity in glassy systems fl5l - l2lj ]. Therefore, a theoretical 
method for describing fluctuations near a saddle-node 
bifurcation might be useful for studying wider systems 
including glassy systems (see Sec. [V] as a related discus- 
sion). 

With this background, in the present paper, we ana- 
lyze a Langevin equation for a quantity (j) in which the 
deterministic part undergoes a saddle-node bifurcation. 
Especially, we investigate statistical properties of relax- 
ation behavior near the bifurcation point with small noise 
intensity T. Note that T is proportional to the inverse of 
the number of elements in cases of many-body systems 
with an infinite range interaction Q or defined on a ran- 
dom graph [TTI. [T3] . Then, the fluctuation intensity x<p(t) 
of relaxation trajectories exhibits a divergent behavior 
similar to those observed in the dynamical heterogene- 
ity. The aim of this paper is to describe this divergent 
behavior theoretically. 

The main idea in our theoretical analysis is to express 
a single trajectory in terms of a singular part and the 
others. Concretely, we focus on the exiting time 9 from a 
saddle-point and find that 9 exhibits the divergent behav- 
ior ({9 - (9)f) ~ T~ 2 / 3 in the limit T -> with e = 
fixed and ((9 - (9)) 2 ) ~ e~ 5 / 2 in the limit e -> with 
T« 1 fixed. In addition to these divergent behaviors, 
we derive a statistical distribution of 9, by which x^(£) is 
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calculated. This idea does not only make the calculation 
possible, but also provides us an insight for the nature of 
the divergent fluctuations near a saddle-node bifurcation. 
That is, the most important quantity that characterizes 
the fluctuations near a saddle-node bifurcation is the ex- 
iting time from the saddle-point. 

This paper is organized as follows. In Sec. UH we 
present a model, and display numerical results on di- 
vergent behavior of X<^(*)- I n Sec. Mil on the basis of a 
simple phenomenological argument, we derive critical ex- 
ponents characterizing the divergent behavior. We also 
present the basic idea of our theory. In Sec. IIV1 by em- 
ploying a singular perturbation method, we construct a 
systematic perturbation theory so as to determine the 
statistical properties of the important quantity, 9. Sec- 
tion [V] is devoted to concluding remarks. 



II. MODEL 

Let <f)(t) be a time dependent one-component quantity. 
We study relaxation behavior described by a Langevin 
equation 

9 t = / e (0)+£. (1) 

Here, / e (0) is assumed to be 

/ e (0) = -e0-0(0-l) 2 , (2) 

where e is a small parameter. £ in (TTJ) represents Gaussian 
white noise that satisfies 

(£(t)£(0>=2T*(t-f), (3) 

where T represents the noise intensity. A qualitative be- 
havior of deterministic trajectories for the system with 
T = is understood from the form of a potential defined 
by f £ = —d^v e , where the potential v e ((f>) is written as 

« e M = §0 a + 5* 2 -§* 8 + ^ 4 . (4) 

As displayed in Fig. [1] there is a unique stable fixed point 
= when e > 0, while there are two fixed points = 
and — 1 when e = 0. Since the latter is marginally 
stable, it is called a marginal saddle. Furthermore, when 
e < 0, there are two stable fixed point and one unsta- 
ble fixed point. This qualitative change of deterministic 
trajectories at e = is called saddle-node bifurcation. 
Throughout this paper, we assume the initial condition 
(f)(0) = 0o > 1 and e > so that the trajectories pass the 
marginal saddle at = 1. 

We denote a trajectory 0(f), < t < oo, by [0]. All 
the statistical quantities of trajectories are described by 
the path probability measure 

^4> 
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FIG. 1: (color online). Potential v e (cj>) for a few values of e. 

where is the normalization factor which is indepen- 
dent of [0], and .F([0]) is determined as 

H[4>]) = lf dt ~ fMf + 2T/ e '(0)) . (6) 

(See Appendix [Al for its derivation.) The last term in ([6|) 
corresponds to the Jacobian associated with the trans- 
formation from a noise sequence [£] to the corresponding 
trajectory [(f)]. 



A. Numerical simulations 

Before entering the theoretical analysis of ([T]), we re- 
port results of numerical simulations. The Langevin 
equation was solved numerically with the Heun method 
[22J with a time step 10 ~ 4 , and the initial condition was 
fixed as 0(0) = 0o = 1.2. The expectation value (A) of 
a fluctuating quantity A was estimated as the average of 
one-hundred data. By using ten independent samples of 
the estimated values, the value of (A) is conjectured with 
error-bars. 

In Fig. [2 we display nine trajectories (f>(t) for the sys- 
tem with e = and T = 10~ 6 , where each trajectory is 
generated by a different noise sequence. The trajectories 
are clearly distinguished despite a rather small value of 
T . It is seen that a major difference among the trajecto- 
ries is the exiting time from a region near = 1. Similar 
behaviors are observed for other small values of e and T. 

In order to clarify (e, T) dependence of the relaxation 
behavior, we investigate (4>(t)) for several values of e and 
T. As is seen from Fig. G2 (0(f)) exhibits the two steps 
relaxation, and the plateau regime around ~ 1 becomes 
longer as T is decreased with e = fixed (see Fig. [3]), or as 
e is decreased with T = 2~ 20 fixed (see the inset of Fig. [3]). 
These qualitative behaviors are easily conjectured from 
the form of the potential v e ((f>). 

An impressive feature of the relaxation behavior is 
qualified by the fluctuation intensity of 4>{t): 

X4> (t)^T-\(^(t))-(cp(t))\ (7) 
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Note that X<p(t) 1S independent of T in the limit T — >• 
when fluctuations are not singular. As displayed in Fig.[H 
each X$(t) for small e and T takes a maximum value at 
a time i*. Furthermore, it is seen from these graphs that 
both the time i* and the amplitude X4>{t*) increase as 
T is decreased with e = fixed (see Fig. 0]) or as e is 
decreased with T = 2~ 20 fixed (see the inset of Fig. 
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FIG. 2: (color online). Trajectories <f>(t) with e = and 
T = 2 -20 . The nine trajectories are generated by different 
noise sequences, respectively. 
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FIG. 4: (color online). X<#>M f° r several values of T. 
Inset: X </>(£) f° r several values of e. T — 2~ 20 . 



X<t>(t*)T seem to coincide with each other in the regime 
£ y-2/3 ^ ^ for different small values of T. The graphs 
in the insets of Figs. [5] and [6] also indicate the relations 
([ID]) and (HH) in the limit T — > with small e fixed. 
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FIG. 5: (color online). t,T 1/3 as functions of eT~ 2/3 for sev- 
eral values of T. Inset: t* as functions of e. The guide line 
represents a power-law function t* ~ e~ 
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FIG. 3: (color online). (<j>(t)) for several values of T. e = 0. 
Inset: (4>(t)} for several values of e. T = 2~ 20 . 

Now, we describe the divergent behaviors of t* and 
X(f>(t*) quantitatively. We first note that the behaviors 
are classified into two regimes, (i) e < T 2 / 3 < 1 and (ii) 
T 2/3 < e < 1. We observe 

U * T' 1 / 3 , (8) 
X4,(U) * T~\ (9) 

in the regime (i), while 

t. - e" 1/2 , (10) 
X,(t.) ^ e~ 5/2 , (11) 

in the regime (ii). The relations © and ([9]) are con- 
jectured from the graphs in the mainframes of Figs. [5] 
and [51 respectively. Indeed, all the data of i^T 1 / 3 and 




FIG. 6: (color online). )T as functions of eT" 2/3 for 

several values of T. Inset: x<t>(t*) as functions of e. The guide 
line represents a power-law function x<p(t*) — £~ 5 ^ 2 - 

The specific purpose of this paper is to provide a the- 
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oretical understanding of the divergent behaviors given 
by (JgJ), ©, (UHl), and $TB). Here, we note that (cf>(t)) ex- 
hibits the discontinuous behavior from e = to e — > 0+ 
with T = fixed or from T = to T -> 0+ with e = 
fixed. The coexistence of the discontinuous nature of 
(4>(t)) and the critical nature of x<£(^*) i s a character- 
istic feature of stochastic dynamics near a saddle-node 
bifurcation. Such a coexistence, which is called a mixed 
order transition, has been observed in several systems in 
the context of dynamical heterogeneity (15h21|. See a 
related discussion in Sec. |V] 
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III. PHENOMENOLOGICAL ANALYSIS 

In this section, we present a phenomenological argu- 
ment for deriving the relations ©, ©, (fTO]). and ([TTj). 
This argument also provides an essential idea behind a 
systematic perturbation method which will be presented 
in Sec. |lVl 

We first note that there are two small parameters e and 
T in this problem. Despite of this fact, the discontinu- 
ous nature of (</>(£)) causes difficulties in the theoretical 
analysis. Indeed, we need a special idea so as to develop 
a perturbation method. Let us recall the typical trajec- 
tories displayed in Fig. [5J All the trajectories are kinked 
in the time direction. Here, the kink position, which cor- 
responds to the exiting time from a region around <j) — 1, 
fluctuates more largely than the other parts of the tra- 
jectories. These observations lead to a natural idea that 
a kink-like trajectory is first identified as an unperturbed 
state. 

We illustrate the idea more concretely. Let us consider 
two special solutions of dt4> = fo(4>)- The first is the 
special solution ^> sp (t) under the conditions </> sp (t) — > 1 as 
t — > — oo, <p S p{t) — > as t — > oo, and (f> sp (0) = 1/2. 4> sp (t) 
takes a kink-like form that connects the two fixed points. 
(See Fig. [71) The second one is the special solution <pBo(i) 
under the conditions </>bo(0) = <fro and 4>Bo(t) — > 1 as £ — > 
oo. (t>Bo(t) connects the initial value and the marginal 
saddle. (See Fig.[7J) Then, by introducing a time 9 which 
corresponds to a kink position, we express trajectories as 

d>(t) = sp (t -9) + (0bo(<) - 1) + <p(t - 0), (12) 

where <p(t — 9) represents deviation from the superposi- 
tion of the two special solutions. 

Now, fluctuations of <f>{t) are expressed in terms of 
those of 9 and <f. Then, it is reasonable to assume that 
if does not contribute largely to the divergent part of 
the fluctuation intensity of 4>{t). With this assumption, 
{4>(t)) and x<j>(t) are estimated as 

(0(<)>=<M*-0)>. (13) 

and 

x*(*) - ^ (<M* - of) - fa*(t - o)f) - (1 4 ) 



FIG. 7: (color online). Functional forms of sp (i — 9) and 
<f>Bo(t)- 

where the statistical average is taken over 9 with a distri- 
bution function of 9. On the basis of these expressions, 
we provide a phenomenological argument by which we 
can determine the exponents characterizing the divergent 
behavior of t* and x<f>(t*)- 

The argument is divided into three parts. In Sec. MI Al 
we first derive the scaling forms of (9) and ((9 — (9)) 2 ) 
in the two regimes, T 2/3 < e < 1 and e < T 2/3 < 1, 
respectively. Second, based on these expressions, we 
conjecture the distribution functions P(9) for the two 
regimes. Finally, in Sec. IIII Bl by using P{9), we calcu- 
late the critical exponents of and X</>(i*), which were 
observed numerically in ©, ©, (ITU1) . and (fTT]). 



A. Statistical properties of 6 

The fluctuation intensity of 9 is defined as 

Xo = ^{{9 2 ) {9)\ (15) 
We then assume scaling relations 

(9) = T-O'f^T-^'e), (16) 
Xe = T-t'/'-MT-V'-e), (17) 

for small e and T, with new exponents and 
7'. Here, the scaling functions fi(x) and f2{x) satisfy 
two conditions: (i) /i(0) and /2(0) are finite and (ii) 
fi(x) ~ and /a(^) — x" 1 for x ^> 1. The latter 
condition implies that (9) and xe are independent of T 
in the regime T 1 ^* -C e <C 1, where 9 is expected to 
obey a Gaussian distribution with variance proportional 
to T (see (27])). 

We shall determine the exponents , v*, and 7' in (fT6|) 
and (|17[) . Since a characteristic time around the marginal 
saddle is expected to obey the same scaling relation as 9, 
we investigate the local behavior near <f> — 1. Concretely, 
we substitute <f> = 1 + <f into ^ and ignore higher powers 
of f from an assumption that \tp\ is small. We then obtain 
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the probability measure for [ip] as 

P(M)=exp (-HM)/T), 



(18) 



where 



1 



HM) = 4j dt [(( dt( P + e + t P 2 ) ~ 4T( P) + const{19) 

From this expression, when T = 0, we immediately find 
that ip has a scaling form ip(t) = e 1 ^ 2 <pi(e 1 ^ 2 t). Since the 
characteristic time scale in this case diverges as e -1 / 2 , we 
find that £' = 1/2. On the other hand, when e = 0, we 
obtain another scaling form y>(t) = T 1 / 3 ^^ 1 ^) which 
leads to C' l v * — 1/3- We thus derive — 3/2. From 
this result, we expect that a distribution function of is 
expressed as a T- independent function of 0T 1 / 3 when e = 
0. This expectation leads to a relation 7'/ 24 = 2Cf / f* + l, 
which yields 7' = 5/2. These results are summarized as 



(0) = T- 1 / 3 / 1 (eT- 2 / 3 ) 
Xe = T-^f 2 (eT-^) 

which, in particular, involve the results 

(9) ~ T- 1 / 3 , 

in the regime e <C T 2 / 3 <C 1, while 

(0) 
Xe 



,-5/2 



(20) 

(21) 

(22) 
(23) 

(24) 
(25) 



in the regime T 2 / 3 -C e -C 1. As shown in Figs. [8] and GO 
numerical results are consistent with (|2"0")) and (|2ip . 
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FIG. 8: (color online). (9) T 1/3 as functions of eT _2/3 . The 
guide line represents (0) = 2 ^6 ly ^ e~ 1//2 , which is our theo- 
retical result (|131[) . The graphs for six different values of T 
tend to converge to one curve as T is decreased. 
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FIG. 9: (color online). xeT 5/3 as functions of eT~ 2/3 . The 
guide line represents \e = 2 ^6 1,/4 ^ e _5//2 , which is our theo- 
retical result (|132[) . The graphs for six different values of T 
are on one curve. 



e = and T <C 1, as the representative in the regime 
e <C T 2 / 3 <C 1. In this case, an exiting event from a re- 
gion near 0=1 occurs by effects of small noise. We then 
naturally expect that trajectories reach to cj> = 1/2 ran- 
domly as a Poisson process except for a short time regime 
in which the behavior depends on initial conditions. That 
is, we conjecture that 9 obeys the Poisson distribution 
when 9 is much larger than some cut-off value 9 cut . Fur- 
thermore, from (|2"2"j) and (|2U)) . the distribution function 
P(9) is expressed as P(9) = P(9T 1 / 3 ). These considera- 
tions lead to 



P(9) = 



rpl/3 e -a9T 1/3 



(26) 



for the regime 9 > cu t- The positive constant a and 
the normalization constant Zg 1 ^ are independent of T. 
We also assume that cut = 0(T -1 / 3 ), which means 
/ e ~ t d9P{9) f 1 in the limit T -> 0. We thus intro- 
duce a T-independent parameter as cut = cut T 1 / 3 . The 
expression (l26t will be derived in Sec. IIV1 

Next, we consider the regime T 2 / 3 < e < 1. When 
e =/= and T = 0, P(0) is the delta function 5(0 - (0)), 
because there are no fluctuations. For sufficiently small 
T, the distribution function is smeared slightly, and this 
leads to a conjecture that obeys a Gaussian distribution 



P{9) = 



1 _ i ce- 
re T ' 



: 



,(2) 



(27) 



(2) 

where Z s is a normalization constant. The expression 
([27]) will be derived in Sec. HV1 



Based on these results (g2j, (@3j, and ((25j) , we 

conjecture functional forms of P(0) from which we will 
calculate x^W- 

First, we consider the regime e <C T 2 ' 3 < 1. In or- 
der to simplify the argument, we focus on the case that 



B. Calculation of x^t) 

On the basis of the results (|2l)l) and (|2"T1) , we calculate 
the critical exponents of £» and X<t>{t*)- Here, for conve- 
nience of calculation, we introduce the Fourier transform 
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of <p sp (t) as 



o^l) = I ^ sp ( w )e Mi . 

Z7T 



(28) 



Then, from the assumptions (fT3")) and ([T4")). we estimate 
and 



^0 sp ( W )e- t (e— 9 ), (29) 



<'■■>< '> 2 > = / ^ / f^M^'K^'^ 



-z(k;-t-u/)(? 



By substituting ([2T)[) and (|50|) into (JTJ), we obtain 



1 f du) f du' ~ 



n e - i(u+u ' )e ^ _ L-io/eyj (31) 

In the following paragraphs, we shall calculate x<t>{t) f° r 
the two regimes T 2/3 < e < 1 and e < T 2/3 < 1, 
respectively. 

First, we consider the regime e <C T 2 / 3 <C 1 with set- 
ting e = 0. By using the distribution function ([2^]). we 
have 



d6e- %uje P{6) 

T l/3„-a<9T 1/3 - 

de— 



Al) 



(32) 



We assume that the first term of ([32]) is neglected when 
T-C 1. Then, by performing the integration, we calcu- 
late 



T l/3 



AD 



TV3 (a + zT-VSo;) 



(33) 



Similarly, we obtain 



i(cj+iy)#\ _ 



T l/3 

7« 



-e cut T 1 / 3 (a+i(w+c l ;')T- 1 '' 3 ) 



TVs ( a + iT-V3 (u, + a /)) 



(34) 



Furthermore, since 4> S p(t) shows a quick change from 
</'sp(i) = 1 to 4> S p(t) — around the kink position t = 9 1 
we approximate ^ S p(i) as 0(— t), where 0(i) is the Heav- 
iside step function. With this approximation, the Fourier 
transform of </> sp (t) becomes 



(f) sp (uj) = lim 



-1 



&—*•□+ iuj — b 



(35) 



Here, we introduce a scaled time t = tT 1 / 3 . We then 
substitute flM]), (JMJ), and ([33 into §1$ . The result is 



3 i(w+o;')T- 1/3 (t^e cut ) 



T I 2tt I 2n b->o+ iuj — b iuj 1 



-S cut a 



(a + iT-V3 (w + a;')) 



-20 cut a 



(a + tT-VSw) ( a + iT-V3 w /) 



V 



(36) 



(30) By the transformation of integrable variables, w = 



jT- 1 / 3 and 



'T 1 / 3 , X4>(t) is expressed as a scal- 



ing form x<t>{t) = X<t>(T l/3 t - 9 cut )/T. This implies that 
U = 0(T-^ 3 ) and \<t>(t*) = 0(T- r ) in the limit T —> 0. 
Furthermore, we expect that these exponents are valid 
for e and T that satisfy e <C T 2 / 3 < 1. In this manner, 
we have obtained the results consistent with © and ©. 

We next consider the regime T 2 / 3 < e < 1 by using 
the distribution function ([27]) . The Gaussian integration 
with respect to 9 leads to 



<e-^ e ) 



and 



-i{pj-\-uj )9 



e -i(^W)W e - X9(w+ / )2T . 



(37) 



(38) 



Then, by substituting ([37]) and ([38]) into ([29]) and (p3Tj) . 
we obtain 



(0(0) 



^ sp ( W )e^«)e-^, (39) 



and 



Xfl (^ 2 +u/ 2 )T 



(40) 



From these, we derive 

X<t> 



1 / fe — 1 \ 

W = »Em M ~ 9t & <#)> ■ (41) 
fe=i ' ^ ' 



Here, let t w = VxoT be the width of the distribution 
of 9. When t w -C 1, the expressions ([5^[) and (|4"T[) arc 
further simplified. Indeed, we may estimate 



and 



X*(t) = Xe(d t <f> sp (t - (9))) 2 . 



(42) 



(43) 
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These results show that x<t>{t) takes a maximum at t = 
£*, where i* ~ (0) ~ e^ 1 / 2 . Furthermore, we obtain 
X<t>(t*) — Xo — e~ 5 ^ 2 - These results are consistent with 
(fT0|) and pip . Note that they are valid in the regime 
t w <C 1 which means T 2 / 5 « 6 « 1. It seems that there 
is no power-law behavior in the regime 
In fact, Fig. |5] suggests that X^(**) does not converge to 
one universal curve in the whole region. 

To this point, we have explained that the singular be- 
havior of X4>(t) is determined by the statistical distribu- 
tion of 9. Our analysis shows that 9 is the most important 
quantity for characterization of the divergent fluctuations 
near the saddle-node bifurcation. This claim is also con- 
jectured from a fact that the statistical properties of 9 
are simpler than those of (j>. (Compare Figs. [5] and [6] 
with Figs. [5] and [HI) Thus, we have focused our theoreti- 
cal analysis on the derivation of the statistical properties 
of 9. Note that the scaling relations (|2"2")l and (1241 were 
derived in Refs. [!, [Hj], and arguments closely related 
to ([25)1 and (f2"5)) were also presented in Refs. [all]. in 
particular, all the statistical properties of 9 are described 
by the analysis of the backward Fokker-Planck equation 
23]. However, in the previous approaches, a perturba- 
tive calculation with small e and T seems quite compli- 
cated. In such cases, it would be almost impossible to 
analyze spatially extended systems. In order to improve 
the situation, in the next section, we develop a system- 
atic perturbative calculation for the distribution function 
of 9 within the path-integral formulation ([5]) with ([5]). 



IV. ANALYSIS 

Our theory basically relies on the idea mentioned in 
Sec. IIIII That is, we start with the expression (|T2"j) and 
derive the distribution function of the exiting time 9. 
Formally, the derivation might be done by performing 
the integration of -P([0]) with respect to <p with 9 fixed. 
However, as far as we attempt, it seems difficult to carry 
out this integration by a standard path integral method. 
One difficulty originates from the existence of a tran- 
sient region before passing the marginal saddle. This 
contribution is described by the interaction of 0bo and 
9, and yields a non-trivial distribution of 9 in the regime 
9 <C (9). The other difficulty arises in calculation of a 
perturbative expansion around the solution <p sp . Since 
the solution approaches the marginal saddle in the limit 
t — > — oo, the stability of the solution is marginal. In 
such a case, a naive perturbation induces a singularity. 
We thus need to reformulate the perturbation problem. 

In order to overcome the two difficulties, in this pa- 
per, we employ a method of fictitious stochastic processes 
[24j [. Concretely, we introduce a variable <fr(t,s) with a 
fictitious time s and define a fictitious Langevin equa- 
tion whose s-stationary distribution function is equal to 
the path probability measure given in ([5]) . The Langevin 



equation is written as 



d a (j)(t, s) 



6F([ 



+ V2Tr](t, s), 



with 



( V (t,s)r)(t , ,s , ))=6(t-t')5(s-s'). 



By substituting (|6]) into (|44|) . we write explicitly 



d s 



(44) 



(45) 



(46) 



where 



Fe,T{4>) 



1)(0-1) 3 + J T (<I>) + G e (<t>)(47) 



with 



M4>) 



-(6^-4), 



1)(20-1) 



(48) 



By interpreting t as a fictitious space coordinate, we re- 
gard (|4"6"| as a reaction-diffusion system with the bound- 
ary condition 0(0, s) = 4>q. Particularly, since the system 
is bistable, one may employ techniques treating kinks in 
such systems [H, [26[. As a result, as will be shown in 
subsequent sections, the interaction between 0bo and 9 
can be formulated as a perturbation and the problem 
arising from the marginal stability of sp can be treated 
in a proper manner. We note that interesting noise ef- 
fects in kink dynamics in bistable systems were reported 
in Ref. [27|. 

As discussed in Sec. lU the qualitatively different be- 
haviors were observed depending on the regime either 
e < T 2 / 3 « 1 or T 2 / 3 < e < 1. Correspondingly, the 
perturbation theory is developed for each regime. Since 
the basic idea behind calculation details is in common 
to both the regimes, we provide the full account of the 
perturbation theory for the regime e <C T 2 / 3 < 1 in 
Sec. HVAl ITVBl and ITVCl Then, we discuss briefly a 
perturbation theory for the regime T 2 / 3 < e < 1 with 
pointing out the difference from the regime e <C T 2 / 3 <C 1 
in Sec. ITVDl 



A. Formulation 

1. Unperturbed system 

Since we focus on the regime e -C T 2 / 3 <t^. 1, one may 
choose the system with e = and T = as an unper- 
turbed system. However, we cannot develop a perturba- 
tion theory with this choice. In order to explain the rea- 
son more explicitly, we define a potential function Vr{4>) 
by 



=0,T 



(0) = -dfiVrtt), 



(49) 
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FIG. 10: (color online). Potential Vt(4>) for T = and T = 
0.01. 
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FIG. 11: (color online). Unperturbative potential V u ((j>), and 
V T (4>) for T = 0.01. 



where Vr(l) = 0. This potential is calculated as 

Vrtf) = --^f^ ~ !) 4 + ~ 1)(3^ - !)• (50) 

In Fig. [TOl the functional forms of Vt (</>) are displayed 
for T = and T = 0.01. Here, the special solution 
fa p (t) in (fT2"j) connects the two maxima of Vr=o (</>)■ How- 
ever, since the curvature of the potential Vr=o( ( / > ) at 
the maximal point <f> = 1 is zero, a perturbative cor- 
rection to the solution fa p (t) exhibits a divergent behav- 
ior. (See an argument below (|88p.) In order to avoid 
the singularity, we must choose an unperturbed poten- 
tial different from Vt=q (</>)• Based on a fact that Vt(4>) 
has the two maxima at <p = fa = —AT + 0(T 2 ) and 
(j) = fa = 1 + T 1 / 3 + 0(T 2 / 3 ), where V T {fa) = T/2 and 
V T (fa) = 3T 4 / 3 /4 + 0(T 5 / 3 ), we define an unperturbed 
potential V n by the decomposition 



Vrtt) = V u (fa+A(fa, 



(51) 



where A(c/>) is chosen such that the two maximal values 
at <p = fa and <p — 4>2 of Ki(0) are identical; that is, 
V u (fa) = K(^) - T/2 and V u '(^) = V^fc) = 0, as 
shown in Fig. 1111 Furthermore, in order to have a simple 
argument, we impose a condition that the curvature at 
each maximum of V u is equal to that of the corresponding 
maximum of Vt with ignoring contribution of 0(T 4 / 3 ). 
These conditions are satisfied by setting A(<^i) = 0, 
A(&) = M<h) ~ V T (<t>i) = 3T 4 / 3 /4 - T/2 + 0(T 5 / 3 ), 
A'(fa) = A' (fa) = 0, A"{fa) = 0(T 4 / 3 ), and A"(fa) = 
0(T 5 ' 3 ). As an example, one may choose 



A(0) 



3 T 4/3 _ ^ 

4 2 



0(T 5 / 3 ) (12 Ue=o (0) + O(T))(52) 



By using the decomposition (|5 1 [) and defining T u = 
— ^Ki: w e rewrite (|4"o]) as 

9,0 = ia t 2 0-F u (0)-G e (^)+A'(0) + V2T77.(53) 

In our formulation, we treat the last three terms as per- 
turbations. 



2. Expression of solutions 

Since we have replaced the unperturbative potential, 
we reconsider an expression of trajectories. We first de- 
fine the special solution fa, of the unperturbed equation 



Id 2 



2 dt 2 



F u (fa) =0 



(54) 



under the conditions fa(t) — > fa when t — > — oo and 
fa(t) — > fa when t — > oo. We also impose fa(0) = 1/2 
in order to determine fa uniquely. This solution corre- 
sponds to the kink solution in the real time direction and 
describes the relaxation behavior from <j> = fa(~ 1) to 
(j) = fa(pi 0). The functional form of fa can be obtained 
by the integration of 



d t fa = -2^V u (fa) - V u (fa). (55) 
In Fig.[H we show fa(t) for T = 10~ 3 and 10~ 5 . 




FIG. 12: (color online). </>* as functions of t. 
Next, let fa be the solution of 







(56) 



2 <% 2 

under the conditions that fa {t) — > fa as t — > oo and that 
^b(0) = fa- This solution describes a typical behavior of 



9 



(f> near t = when T is sufficiently small. By using these 
two solutions, we express the solution of ([5^)1 for a given 
r\ as 

<f>(t, s) = Mt~ 6(a)) + (Mt) ~ <h) + pit - 0(a), s)(57) 

where p represents a possibly small deviation from the 
superposition of the two solutions. 



3. Linear stability analysis 



3>n 




0.6 



0.4 



0.2 



-0.2 



1 „ 1 1 
T= 10"^ 


i 


i 


T= 10" b 






F u '(**(u)) 






i i i 


i i 



FIG. 14: (color online). $0 as functions of u. 
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FIG. 13: (color online). F^((j)*(u)) as functions of u. 

As a preliminary for a systematic perturbation theory, 
we perform the linear stability analysis of Hereafter, 
we set u = t — 9. The stability of the solution </>* is 
determined by eigenvalues of the linear operator L given 
by 



2du 2 



L /•>,(«))• 

Let us consider the eigenvalue problem 
L<S>{u) = -A$(u). 



(58) 



(59) 



This problem is equivalent to an energy eigenvalue 
problem in one-dimensional quantum mechanics, where 
Fu(cj>*(u)) and A correspond to the potential and an en- 
ergy eigenvalue, respectively. The graphs of i 7 ^ (</>*(«,)) 
are shown in Fig. 1131 The asymptotic behaviors are cal- 
culated as F' u (<j)*(u)) -> 1/2 + 0(T) as u -> 00, and 
K{<P*{ U )) -> 3T 2 / 3 + O(T) as u -> -00. 

Since 0*(i — 0) for any is a solution of (|5"4"]) . we may 
take the derivative of ([54| with respect to 9. We then 
have 



0. 



(60) 



This implies that there exists the zero-eigenvalue, for 
which the normalized eigenfunction $o(u) is determined 
as 



(61) 



where 



du(du4>* ( U ) Y 



(62) 



Here, by using d u 4>* — f e —o((f)*) + 0(T), we obtain 



r = 



#/ e=o (0) + O(T) 



v e = mi + o(T) 



= 12 + ° {T) - 



(63) 



The zero-eigenfunction &o(u) corresponds to the Gold- 
stone mode associated with a time translational symme- 
try. The functional form of &o(u) is shown in Fig. [14] 
Since there is no node in the u profile of $o(u), the 
minimum eigenvalue must be zero as in quantum me- 
chanics. Since the height of the potential in the limit 
u — > —00 is 0(T 2 / 3 ), it is expected that there is no 
other discrete eigenvalue when T is sufficiently small. 
Next, we consider continuous eigenvalues for A > A m = 
lim^ooF^O^u)) = 3T 2 / 3 + 0(T). The correspond- 
ing eigenfunctions are characterized by the asymptotic 
plane waves with the eigenvalues A = A m + fc 2 /2 and 
A = (1 + fc a )/2, where /cr and &l are non-zero real num- 
bers representing wavenumbers of the asymptotic plane 
waves, in the limit u — > —00 and u — > +00 respectively. 

We denote by <&\{u) the eigenfunction corresponding 
to the eigenvalue A. Since the eigenvalues greater than 
1/2 are degenerated, $\(u) ^ $\(u) for A > 1/2, where 
* represents the complex conjugate. Here, it is conve- 
nient to introduce $_x(u) as <&-\(u) = §>\(u). Re- 
lated to this introduction, we also define the index set 
A = {A|A < —1/2, A > A m }. Then, we may choose the 
set of eigenfunctions so as to satisfy the orthogonality 
condition 



and 



* («) - (c\A)/vr, 



du$X,(tt)* A (iO = *(A-A'), 



du$ (u)$ A (u) = 0, 



(64) 



(65) 
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for any A, A' e A. Furthermore, we expect the complete- 
ness condition 

*o(«)*o(«') + / d\$l(u)$ x (u') = 6(u-u'). (66) 

JA 

By using these eigenfunctions, we expand p[u, s) as 



p(u,s) = / d\iP x (s)<P x (u), 

J A 



(67) 



where the contribution of the zero eigenfunction is not 
taken into account in the expression (|67j) so that the ex- 
pression (|57p is uniquely determined. 



B. Perturbation theory 



By substituting (|57| into (|53|). we obtain 



-fld u (j)* + d s p 

-G e (^) + A'(^) 
+L(0„)(M*) - <fe) + V2T?y + 



(68) 



i. Lowest order result 

We start with the calculation of f^ 1 / 2 ) and p^ 1 / 2 ). We 
substitute the expansions (|70|) and (|7T|) into f[68]) , arrange 
terms according to powers of p, and pick up all the terms 
proportional to 0(p 1 ^ 2 ). We then obtain 

- Q^d u ^ + = L(^)p(V2) + ^rj, (72) 

We rewrite (fT^j) as a linear equation for p^ 1 / 2 ) with the 
expression 

(d s - Lifa)) =B(4>*, V ). (73) 

Because L(<^*) possesses the zero-eigenvalue, it is not the 
case that there exists a unique bounded p^ 1 / 2 ). That is, 
there exists no bounded solution or there are an infinitely 
number of solutions. In order to proceed the calculation 
further by obtaining p^ 1 ' 2 ' , we impose the solvability con- 
dition under which the latter case is chosen. The solv- 
ability condition in this case is written as 



where we set 



This yields 



where 



du® (u)B((p*,r)) = 0. 



rf2 (1/2) = \/27T77e, 



(74) 



(75) 



(69) 



We consider a perturbation expansion of p and 57 with 
respect to a small parameter. In the present problem, 
G e (</>*), A'(0»), L(<f>*)(<i)B(t)-<h), an d V^Trj are treated 
as perturbations. In order to formulate the problem con- 
cretely, we introduce a formal expansion parameter p in 
front of <3 £ (0»), A' (<£»), L(0,)(0b(*) - 2 ), and T in the 
noise term. We solve the equation perturbatically 
under the assumption that p and fl can be expanded in 
p. Note that p=fi=0 when p = 0. Concretely, due to 
the small noise term which is 0(p 1 ^ 2 ), we assume 



(70) 
(71) 



As we will see below, all the terms of p^ 1 ' and fiW can 
be calculated in principle. Here, it should be noted that 
p is not directly related to small parameters e and T. 
Therefore, for example, if one wishes to derive the solu- 
tion valid up to 0(T 2 ), we do not have a quick answer 
to the question how many orders of expansion in p are 
necessary. Although this aspect makes the analysis com- 
plicated, we will find that the formulation leads to a sys- 
tematic expansion. These preliminaries presented in this 
section are standard in a singular perturbation method 
[Hill]. With this setting up, we calculate p w and 
in sequence. 



du(d u 4>*)r){u + 9,s)/Vr. (76) 



We note that ([76]) satisfies 

( m (s)r,g(s'))=S(s-s') 



(77) 



Furthermore, under the solvability condition, one can 
determine statistical properties of p^ 1 / 2 - 1 from (f73|) . In a 
manner similar to (|67j) . we expand p*- 1 / 2 ) in terms of the 
eigenfunctions of L(<j)*). Then, for any non-zero eigen- 



value A, the coefficient i^^ 2 ^ 



with 



obeys a Langevin equation 
(78) 



-A^ 1/2) + VWr lx , 



du$\(u)ri(u + 9, s). 



From this, we obtain ( ipi 1 



= and 
T 



K /2 ^ /2) *) e = ^(A-n 



(79) 



(80) 



where () g denotes an expectation value with 9 fixed. The 
fluctuation intensity of p*- 1 / 2 - 1 is then calculated as 



(P 



(1/2H2 



d\\\<S>x{u)\ 2 

A * 

-TG(u, u), 



(81) 
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where G[z, y) is the Green function defined as 



G(z,y) 



d\-$ x (zM(y). 



Here, we note that G(z,y) satisfies 
1 



(82) 



d 2 z - KiM*)) ) G(z, y) = S(z -y)- 4> (z)$ (y). 

(83) 

Let us estimate G(u, u) in the limit u — > ±oo by defining 

(84) 



m±= lim y/2F^Mu))- 



From the expression of F u (</>*(«)) determined by (|47|) . 
([STll . and ([52l . we calculate 



m+ = 1 + 0(T), 

m_ = y/6T 1/3 + 0(T 2/3 ), 



(85) 
(86) 



in the limit T — > 0. We then define a Green function 
G±(z,y) by 



This Green function is written as 

1 



G ± (z,y) 



m±\z— y\ 



m± 



(87) 



(88) 



We conjecture that G(u,u) approaches G±(u,u) as u — > 
±oo. Then, from (JEE]), we find that ((p (1/2) ) 2 ) e ^ 
0(T 2 / 3 ) for u ->■ -oo, and ((p (1/2) ) 2 ) e ~ O(T) for 
w — > oo. 

Here, we address one remark. If 4> sp (t) were chosen 
as the unperturbative solution, m_ would become zero, 
and therefore p^^ 2 ^(u ~ — oo) would exhibit unbounded 
Brownian motion as a function of s. This singularity 
originates from the marginal stability of sp . This is 
why we choose 0» instead of 4> sp as the unperturbative 
solution, which was mentioned in Sec. IIV A II 



The solvability condition for the linear equation for p( l ~> 
yields 

-TCI* 1 * = + * 2 + * 3 + * 4 , (90) 



where 

*i = — 2 \ du(d u ^)F^)(p(^f, (91) 

/OO 
du(d u ct>*)A'(ct>*), (92) 
-OO 

/OO 
du(d u (f>*)L((f>*)((f> B (u + 6)-<P2), (93) 
-OO 

/OO 
du(d u <t>*)G e (<P*). (94) 
-oo 

^2 is immediately obtained as 



<J/„ 



| - ^T 4 / 3 + 0(T 5 / 3 ). 



(95) 



^4 is also calculated as 



*4 = 



e 2 I^ 2 



+ (0-l) ; 



ieT 2/3 . 



(96) 



Because the calculation steps for and ^3 are much 
longer than "I^ and ^4, we shall present them in the 
subsequent sections. 



3. Calculation of \ti 



2. Next order calculation 

In the lowest order description, the variable 6 exhibits 
unbounded Brownian motion, and therefore it indicates 
a singular behavior. Now, in order to determine the ex- 
ponents characterizing the singularity, we proceed to the 
next order calculation. We substitute (|70)l and (|7T|) into 
and pick up all the terms proportional to O(p). We 
then obtain 

- n^d u ^ + djv = £(^)p (1) - I F ^) (p (i/2))2 

-G e (4>.) + A'{4>.) 
+L(<f>*)(4, B (t)-fa)- (89) 



We calculate <Fi defined in fljJTJ. First, (p (1/2) ) in 

the right-hand side of (|9T]) is replaced with ^(p' 1 ^ 2 '') 2 ^ 7 

because /Z 1 / 2 ) is determined by the linear Langevin equa- 
tion (|78l) . By performing the partial integration and us- 
ing (|58| . the result is expressed as 




+ ^ 2 ((p(V 2 )) 2 ) 

u— — 00 

-\ fjuLd u ((pWf) e . (97) 
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The third term of the right-hand side of (|97|) is further 
rewritten as 



1 



dX- 



x (<f> x (u)d u <I> x (u)* + <P* x {u)d u <l> x (u)) 



i / eL\— 
2 7a A 

1 



-X\$ x (u)\ 2 + -\d u <l> x (u)\ 2 



> a (u)^*a(«)*+*3;(«)^*a(«)) 



(98) 



li— — oo 



By substituting into the third term of (j9"7) and using 
the eigenvalue equation again, we express in terms of 
the difference of boundary values of a quantity H . That 
is, we write 



*! = ff(oo) - H(-oo), 
where H{u) is defined as 



(99) 



If T 
H { u) = --jdX- 



.(100) 



Next, we calculate i?(±oo). In terms of the Green 
function given by (|82|) . we express (|100[) as 

= j [d 2 z G{z,y)-d z d y G{z,y)\\ z=y=u .{lQl) 
Here, it is reasonable to assume that 



Note that r becomes larger as the typical value of 9 is 
larger. Indeed, r — > oo in the limit T — >• 0. The inte- 
gral region [—00,00] in (|93|) should be read as a formal 
expression for [— r, 00] with the limit T — > 0. 

Now, by performing the partial integration and by us- 
ing the relation (|60|) . we obtain 



* 3 = ^d u ^(u) ■ d u (<p B {u + 9) ~ h) 



U — — T 

U — OO 



Since lim^ 
we write 



with 



-z(%M»)) (Mu + e)-4> 2 ) .(105) 

U— — T 

,<9 M (0 B -02) = and lim„_ i . oo (0B-02) = 0, 



* 3 = --(K 1 -K 2 ), (106) 



Ky = (d u <t>*)d u (<t> B -4> 2 )\ u= _ T , (107) 
K 2 = [duidMifo - cj> 2 )\ u= _ T . (108) 

Let us evaluate K\ and K 2 . First, since the T de- 
pendence in this contribution is not singular, we assume 
T — > in this evaluation. Then, 0b and 0* satisfy 
dt<t> = —0(0 — I) 2 - We then derive the asymptotic form 



62 + i 



(109) 



as t — > 00, and 



ff(±°°) = j [d 2 z G±(z,y) - d z d y G±(z,y)] | 2=y=±oc ■ 



(102) 



We then obtain 



*i = H(oo) - H(-oo) 
T 

= — (-2m+ + 2m_) 



0*(m) ~ 02 + 



(110) 



as u 
late 



-00. By using these asymptotic forms, we calcu- 



T 
2" 



(l - V6T 1/3 + 0(T 2/3 )) . (103) We thus obtain 



4- Calculation of $ 3 



^2 



16 

04' 



32 



*3 = TT" 



24 
0*' 



(111) 
(112) 

(113) 



We calculate ^3 defined in (|93|) . First, note that the in- 
tegral region in (|93[) is written in a formal manner. More 
precisely, since the u-integration should be defined in the 
bulk region of the special solution 0*(it), the integral re- 
gion is replaced with [— r, 00], where — r corresponds to 
a matching point between the solutions 0b (u + 9) and 
4>*(u). Since the behaviors of 0b (it + 9) and 0*(it) are 
symmetric around ~ 1, we assume that the matching 
point is 9/2 when 6* ^> 1. From this consideration, we set 

8 , , 

r = -. (104) 



5. Result of Qi 

By substituting (|95j) . (|96|) . (fT03|) . and ([IT5jl into (|90]>, 
we obtain the result of fli as 

We define a potential U(9; T, e) so that (|1 14|) is expressed 
by 

r^i = -d e U, (115) 
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where the potential is determined as 

"2/3/, , 3 



U(6;T,e) 



4 



-eT' 
2 



4 



(116) 



It is worthwhile noting that U(9; T, e) satisfies the scaling 
relation 



U(9; T, e) = TU(9T 1/3 ; 1, eT~ 2/3 ). 



(117) 



We_then define = 9T 1 / 3 , I =_ eT" 2 / 3 , and t/(0;e) = 
[7(0; 1, e). In Fig. El we plot U(9; e) as functions of 
for a few values of e. Here, the first, second, and third 
terms of (II 1 6|) represent the driving force in the negative 
direction of 0, while the last term of (|116[) is the repulsion 
from the boundary = 0. The most probable value of 
0, which corresponds to the minimum of the potential, is 
determined by the balance of these two effects. 




FIG. 15: (color online). U as functions of 9 for a few values 
of e. 



C. Distribution function of 



We combine (|75p and (| 1 14|) . After setting fi = 1, we 
obtain 

Td s 6 = -e 2 \ - \eT 2 ' 3 2 -^_A T ^ + ^ + V2TT Ve , 

(118) 

where r\e satisfies ([77| . Using the potential (|116[) , we 
derive the s-stationary distribution function of as 



P(0; T, e) = — exp I 



(119) 



where Zg is the normalization constant. It should be 
noted that the distribution function satisfies the scaling 
relation 



P(0; T, e) = P(9T 1/3 ;eT- 2/3 ). 



(120) 



The expression (|119p with (111 6[) is the main result of our 
perturbative calculation. 



Unfortunately, the distribution function (jl 19[) with 
(|116p is not the precise expression even in the limit 
e < T 2 / 3 < 1. There is a subtle reason. Recall that 
(|118p is valid up to 0(/i) in the formal expansion se- 
ries. When we calculate higher order contributions in 
the equation for 0, the coefficients in (|118l) are modified. 
For instance, — Tf^ includes the contribution 



(121) 



By estimating this quantity, we have found that this term 
is 0(T 4 / 3 ) + 0(e 2 ) + 0(eT 2 / 3 ). (See Appendix E) The 
determination of the coefficients of these terms of / is 
impossible without the numerical integration. There- 
fore, we did not derive the precise expression in the limit 
e <C T 2 / 3 <c 1. Nevertheless, we present two positive re- 
marks. First, terms of 0(T 4 / 3 ) + 0(e 2 ) + 0{eT 2 / 3 ) do not 
appear beyond some order of expansion in /i. Therefore, 
in principle, one may have the formula determining the 
coefficients of these terms. Second, the scaling relation 
(|120p in the limit e < T 2/3 < 1 seems valid up to all 
orders of expansion in [i. 

In order to check the validity of the scaling relation 
(|120p . we investigate the distribution function P(0;T, e) 
for the case e = by numerical simulations of ([1]). In 
Fig.rjll we plotted P(0; T, 0)T -1 / 3 as functions of 9T 1 / 3 
for several values of T . One may find that four graphs 
for different values of T are not completely collapsed on 
one curve in the region for small 9T 1 / 3 . However, since 
the two curves with T = 10~ 5 and T = 10~ 6 almost 
coincide with each other, we expect that one universal 
curve is obtained when T is decreased further. We thus 
conjecture the scaling relation (|120l) in the limit T — > 
is valid. 

Furthermore, following our theory, we attempt to fit 
these numerical data by assuming the form 



P(0; T, 0) = exp 



1 

T 



iT 4 ' 3 9 - 



b 

03 



(122) 



where c is determined by the normalization condition 
when values of a and b are given. These values well-fitted 
to the numerical data are estimated as ag t = 0.36 and 
tifit = 9, which are compared with our theoretical values 
a* = (2\/6 - 3)/4 = 0.4747 • • ■ and K = 8. The slight 
difference between (afit,^fit) and (a*,6„) comes from the 
contribution of higher order terms than O(fi). (As an 
example, the numerical estimation of I in (|12ip provides 
7~ 0.1T 4 / 3 , which might improve the difference between 
a» and afi t .) 



D. Results for the regime T 2/3 < e < 1 

The theoretical argument developed for the regime e 
T 2/3 < 1 cannot be applied to the regime T 2/3 < e < 1, 
because a deviation p(u) ~ 0(eT^ 1 ^ 3 ) becomes quite 
large there. (See (|B5|) in Appendix |Bj) In particular, 
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FIG. 16: (color online). P{6; T, 0)T _1/3 as functions of 6T 1/A 
for several values of T. The dashed blue line represents the 
functional form (112211 with a = 0.36 and b — 9.0. 



from a fact that the unperturbed potential V n ((j>) goes 
back to the original potential Vt=o in the extreme case 
T = 0, it is obvious that we need to reformulate a per- 
turbation theory. 

The idea is natural and simple: we replace the unper- 
turbative potential V n (cj>) with a potential W n (cj>) which 
is appropriate in this regime. Concretely, instead of (|49|) 
and ([5D]>. we define W e ((f>) by 



F e , T =o(4>) = -d<f,W e ((t>), 
where the functional form of W e (<f>) is given by 

1 ,0 



(123) 



Weft) 



4 

' 4 2 



1) e. (124) 



Here, W e (0) has two maximum points as displayed in 
Fig. 1171 Then, in the manner similar to (|51|) , we consider 




FIG. 17: (color online). W e as functions of (f>. They possess 
two maxima at = and tj> = 1 — e/2. 



the decomposition 

W- e (0) = W u (0)+o(^), 



(125) 



where a((f>) is chosen such that the two maximal values 
of the potential W u ((f>) are identical. With this setting 



up, we study the equation 
1 



(126) 



After that, wc repeat essentially the same procedures 
under the assumption that the last three terms are 
treated as perturbations. Note that the special solutions 
and 4>b are redefined using W u (<fi) instead of V u (<f>). 
Then, in the regime T 2 / 3 < e < 1, we derive the leading 
order expression 



Td s 9 



,1 24 



(127) 



Here, the complicated terms containing p that appeared 
in the analysis in the regime e -C T 2 / 3 <C 1 become 
higher order terms in the regime T 2 / 3 <e<l. Then, 
the s-stationary distribution function of 9 is 



P(0;T,e) = -exp 













f JO) 



(128) 



where Y is the normalization constant. Setting 8 — e 1 / 2 t 
and e = eT~ 2 / 3 , we rewrite P(9; T, e) as P(9; e), where 



Y 



exp 



-3/2 



Since e«l, P(9; e) is further simplified to 



P(9;e) 



1 

7 



exp 



(_^( 



6» - 2(6 



1/4- 



4(6V4) 



(129) 



.(130) 



This is equivalent to (|2"T|) . where (0) and xe are deter- 
mined as 



Xo 



2(f 
2(6 1/4 



l/4 )e -l/2 



-5/2 



(131) 
(132) 



We expect that these expressions are exact in the limit 
T 2 / 3 < e < 1. In Figs. [5] and O we display these 
theoretical results. They are in good agreement with the 
results of numerical simulations. 



V. CONCLUDING REMARKS 

In this paper, we have developed a theoretical frame- 
work for calculating statistical properties of critical fluc- 
tuations near a saddle-node bifurcation. The essential 
idea in our formulation is to choose an unperturbativc 
state in accordance with the bifurcation structure. Since 
trajectories kinked in the time direction have much sta- 
tistical weight near the bifurcation point, we express tra- 
jectories as (j){t) = </>»(t — 9) + (0b(£) — 02) + p(t), where 
4>*(t — 9) represents a one-parameter family of classical 
solutions in the language of the path-integral expression. 
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The parameter 9 is regarded as a Goldstone mode associ- 
ated with the time-translational symmetry. This expres- 
sion naturally provides a divergent behavior, because the 
Goldstone mode is gapless or massless. Indeed, we have 
found that the fluctuation intensity of 9 exhibits the di- 
vergence in the limit e — > with small T fixed and in the 
limit T — > with e = fixed. The divergent behavior of 
X<l>{t) originates from critical fluctuations of 9, and X</>(i) 
becomes complicated due to a non-linear transformation 
to <f> from 9, as shown in Figs. [5] and HO 

Before ending the paper, we wish to explain how the re- 
sults in this paper are related to understandings of other 
apparently different systems. First, our results suggest 
the following general story. When a saddle in a deter- 
ministic description becomes to be connected to an ab- 
sorbing point at some parameter value, fluctuations of 
trajectories exhibit a critical divergence due to the ex- 
istence of a Goldstone mode; and if the saddle is far 
from an absorbing point, the final value of the trajec- 
tory exhibits a discontinuity. Such cases might be re- 
lated to a mixed order transition [l5l - [2ll |. The elemen- 
tary saddle-node bifurcation studied in this paper corre- 
sponds to the simplest one among them. As other types 
of bifurcation associated with a mixed order transition, 
we list up a saddle-connection bifurcation which arises 
in a model for a many-body colloidal system [3(| l3lj |. 
and a mode-coupling transition in a spherical p-spin glass 
model [32l. [33l| . Here, the important message of our paper 
is that one will be able to develop a calculation method 
for the statistical properties of critical fluctuations at a 
mixed order transition by applying the basic idea of our 
formulation to each system under investigation. 

As an interesting and non-trivial, but still simple 
example of mixed order transition, we briefly discuss 
the mode-coupling transition of a spherical p-spin glass 
model with p > 3. The fluctuation property was stud- 
ied by a mode-coupling equation supplemented with an 
external field [34], EH or by the field theoretical anal- 
ysis (3(|. Differently from these previous methods, we 
will be able to consider another theoretical framework 
for critical fluctuations along with the above-mentioned 
strategy. Concretely, we start with a useful expression of 
trajectories as we did. In the thermodynamic limit, the 
equation for the time-correlation function obeys a mode- 
coupling equation [HjEIj], and recently we have derived a 
global expression of the solution near the mode-coupling 
transition |37| . which is similar to (112[) . From this re- 
sult, we may express a fluctuating correlator in terms of 
a Goldstone mode A associated with the dilation sym- 
metry that arises in the slowest time scale. Thus, it is 
natural to describe critical fluctuations near the mode- 
coupling transition in terms of fluctuations of A. The 
concrete calculation will be reported in future. 

The results in this paper also provide some physical 
insights into the dynamics associated with a mixed or- 
der transition even if the calculation has never been per- 
formed yet. For example, let us consider a dynamical 
behavior of a super-cooled liquid (38[- Within a frame- 



work of the mode-coupling theory, a mode-coupling tran- 
sition occurs at some temperature or density [39|. In 
some super-cooled liquids, the behavior associated with 
the transition has been observed approximately [4(J 5l[ . 
The physical picture is well-understood: when the sys- 
tem is near the mode-coupling transition point, a parti- 
cle cannot move freely; and only when the particle gets 
over surrounding particles, it can move. Such an event 
is called an unlocking event [36| or a bond breakage event 
[42l ]. As the mode-coupling transition is approached, the 
frequency of unlocking events becomes smaller and the 
events become correlated more and more in a spatially 
heterogeneous manner, which is called the dynamical het- 
erogeneity. (See Refs. [TH, [l6| as reviews, and Refs. [S- 
|46| as experimental studies, and Refs. 1171 Il8l l42t 

mm 

as numerical observations, and Refs. [19U2ll I34U361 l5J- 
[57| as theoretical studies.) 

From our point of view, we conjecture that an un- 
locking event near the mode-coupling transition might 
correspond to an exiting event from a marginal saddle 
in some local dynamics. Then, the results of our pa- 
per suggest that the most important characterization of 
the dynamical heterogeneity is space-time fluctuations of 
exiting time from a saddle-point. At present, it seems 
difficult to derive such local dynamics theoretically, but 
it is tempting to connect the idea of cooperative arrange- 
ment regions to a marginal saddle in some local dynam- 
ics. Since the distribution function of exiting time from 
cooperative arrangement regions exhibits an interesting 
behavior (58|, the theoretical derivation of the observa- 
tion may be a good starting point for the consideration. 
The final goal in this direction is to obtain a simple ex- 
pression determining space-time fluctuations of the Gold- 
stone mode on the basis of a microscopic particle model. 
We will study it step by step toward this goal. 

Related to the description of dynamical heterogene- 
ity, we remark on a well-recognized conjecture that the 
mode-coupling transition described theoretically is noth- 
ing but a cross-over phenomenon in finite dimensional 
systems. In order to understand the nature of this 
cross-over, one needs to describe an activation process 
from "pseudo" meta-stable states which are not defined 
clearly, but might be connected to that defined in the 
mean-field approximation. Although, physically, the ac- 
tivation process corresponds to a nucleation of some do- 
main, its mathematical expression is highly non-trivial. 
Here, when we apply our analysis to a finite dimensional 
system, in our viewpoint, the cross-over phenomenon is 
equivalent to the finite value of the expectation of the 
Goldstone mode. Thus, we have only to calculate an 
asymptotic tail of the effective potential for the Gold- 
stone mode in the limit t — > 00. Since the analysis of the 
super-cooled liquid is too difficult, we should begin with 
the study of a diffusivel y co upled model of local dynamics 
with ((2]). (See Ref. [l3| as a report on a numerical ex- 
periment with e = 0.) It would be possible to develop the 
mean-field analysis of the spatially extended system, but 
it seems difficult to treat spatial fluctuations accurately 
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even for such a simple system. To develop a systematic 
theory beyond the mean field analysis is a challenging 
problem. 

Finally, let us recall that our formulation is based on 
the fictitious time formalism. One may expect that the 
calculation can be done within a standard Martin-Siggia- 
Rose (MSR) formalism [59] . Such a reformulation is par- 
ticularly important when we study more complicated sys- 
tems. In this context, it is worthwhile noting that the 
third term in (|116p has been obtained in the MSR for- 
mulation with a semi-classical approximation [60(. Such 
calculation techniques in spatially extended systems will 
be developed in future. 
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Appendix A: Path integral expression 

We derive the path-integral expression ([5]) with ([6]) 
from the Langevin equation ([TJ with ([3]). In particu- 
lar, we carefully discuss the derivation of the so-called 
Jacobian term. 

Let At be a sufficiently small time interval. We dis- 
cretize ([T]) as 



,At 



4>n - 4>n-\ = [fe{4>n-l) + fe{<t>n)]-^ 

-K«-iAt + 0((At) 2 ), 



(Al) 



with n = 1, 2, • • • , N. Here, (^n)^^ 1 obeys the Gaussian 
distribution 



/ At \ N/2 
^(fo,£ir-- ,&v-i) = ( ) exp 



At 
4T 



N-l 
n=0 



(A2) 



In the limit At — ¥ 0, <p n is expected to provide <p(nAt) in 
the Langevin equation. 

Let us fix (f>Q. Then, a sequence (£cb£:b ' ' ' j£/v-i) de- 
termines uniquely the sequence ((pi, • • • , 4>n)- Thus, the 
distribution function of the sequence {4> n )n=i i s expressed 
as 



■>N 



) = -Pc(Co,--- ,6r-i) 



9{(f>i, ■ ■ ■ ,4>n) 



(A3) 



Here, the determinant of the Jacobian matrix is calcu- 
lated as 



d{4>i, ■■ ■ ,<I>n 
By using the relation 

JV 



JV 

n 



1 - ~/ £ '(0 n )At 



1 



(At) 



JV ' 



(A4) 



n 



i - \m n )M 



l-Ek'(0n)At + O((At) 2 ) 

n=l L 

ex p(-E^'(0n)At + O((At) 2 )^ 



(A5) 



we obtain 



\4nTAtJ 



N/2 




1 2 



(/c(^n-l) + /e(^))/2 



-2Tf e {4> n )+0{{M)^))) 



(A6) 



By taking the limit At — >• and N — > oo with jVAt fix, 
we write formally ([5]) with (|6"T). 

At the end of this appendix, we remark on a discretiza- 
tion method. One may notice that another discretized 
expression 

K - 4>n-l = fe(<t>n-l)At + At + 0((At) 3 / 2 ) (A7) 

does not yield the Jacobian term, because in this case 
the determinant of the Jacobian matrix 



<9(£o, • • • ,6v-i) 



d{4> 



1 , • ' ' ,<PN 



1 



(At) 



JV 



(A8) 



does not depend on <f>. However, the discretization (|A7I 
provides 



exp 



At 

4T 



N 

E 

n=l 

f>n-l 



1 



(At) N \4nT J 



N/2 



At 



f e (<p n - 1) 



R 



(A9) 



where the term R comes from the product of (j>, 
and the last term (3((At) 3 / 2 ) in the right-hand side of 
(A7) when Q-i ^ s evaluated. Explicitly, R is equal to 
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O((At) 3 / 2 )(0„ - </>„_i)/(At) 2 . We here note that <j> n - 
0„_i ~ ©((At) 1 / 2 ). This leads to R = 0((At)°) = 0(1) 
in the limit At — > 0. Therefore, (| A9|) is not useful in 
the limit At — > 0. With regard to the discretization 
problem, we remark that numerical simulations of the 
discretized form yield ([22]) and (|23]>. too. Here, 

one may confirm that (|2"2"|) and (f2"5)) cannot be obtained 
without the Jacobian term in ([5]) within the analysis of 
the path integral expression. Therefore, this example 
provides an evidence for the claim that the path integral 
expression in the limit At — » always contains the last 
term in ([5]). 



Appendix B: estimation of I 



The first term in (IB2[) is estimated as T£ x 1 /m_ x 
Fq x l/m_, where £ is the length scale over which the 
integral is dominant. By using £ ~ ~ 0(T -1 / 3 ) 

(see ([5o]) h we estimate the first term of (|B2[) as 



T 



dvG(u,v)FZ(Mv))G(v,v) 



~ o(Tx T~ 1/3 x T~ 1/3 x T 1/3 x T -1/3 ) 

~ 0(V/ 3 ). (B3) 

Similarly, the second term in (|B2[) is estimated as 



We estimate the integral I given in (|12ip . Before the 
estimation, we need to evaluate We multiply by 

&\(u) on both sides of and integrate them over 
— oo < u < oo. We then consider the s-stationary state 
and extract the lowest order terms. By using (|67|) . we 
derive 



Thus, we write 



dvG(u,v)G e (0*(w)) 

£ x l/m- x G c (4>^(u)) 

O (V- 1 / 3 x T- 1 / 3 x eT 1 / 3 ) 

O (eT" 1 / 3 ) . (B4) 



d«$3;(tt)^'(^(u)) ((p (1/2) ) 2 ) -°°)> e = °( Tl/3 ) + °( eT ~ 1/3 )- ( B5 ) 



du*J(u)G e (^(«)), 



(Bl) 



where the contribution from the fourth and fifth terms in 
the right hand side of (|89|) have been neglected, because 
they are estimated as higher order terms. From (|67[) . 
(I8T1). and (|82l). we obtain 



T 



dvG{u,v)F^{v))G{v,v) 



+ / dt>G(u,«)G e (0*(t;)). (B2) 



Since the most singular behavior arises around u — > — oo, 
we estimate (p^(u)) in the limit u — > — oo by noting 



From this result, the quantity / in (|121[) is estimated 
as follows. First, (|12ip is rewritten as 

I = -\ J_ l d0»F:(^)(p (1) ("(0*)))^ (B6) 

where u((j>^) is the inverse function of <p*(u). In this in- 
tegration, the most dominant contribution arises from 
the integral region [—(f>2, 1], and this contribution is esti- 
mated as 

7-0 (T 1 / 3 x T 1 / 3 x (T 1 / 3 + eT" 1 / 3 ) 2 ) 

~ 0(T 4 / 3 ) + 0(eT 2 / 3 ) + 0(e 2 ). (B7) 
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